 SUBROUTINE POIS_BOUNDS(nx,ny,nz,dims,nproc,comm3d, &
 & nxpi,nxpe,nypi,nype,nzpi,nzpe, &
 & nxpi_yz,nxpe_yz,nypi_yz,nype_yz,nzpi_yz,nzpe_yz, &
 & nxpi_xz,nxpe_xz,nypi_xz,nype_xz,nzpi_xz,nzpe_xz, &
 & ksize,isize,jsize,nxpi_x,nxpe_x,nypi_y,nype_y)

! use MPI
 use MPI_VARIABLES, only : &
 & decomp

 implicit none

 include 'mpif.h'

 integer, intent(in) :: nx, ny, nz
 integer, intent(in) :: dims(1:3)
 integer, intent(in) :: nproc
 integer, intent(in) :: comm3d


 integer, intent(out) :: nxpi(0:nproc-1), nxpe(0:nproc-1)
 integer, intent(out) :: nypi(0:nproc-1), nype(0:nproc-1)
 integer, intent(out) :: nzpi(0:nproc-1), nzpe(0:nproc-1)
 integer, intent(out) :: nxpi_yz(0:nproc-1), nxpe_yz(0:nproc-1)
 integer, intent(out) :: nypi_yz(0:nproc-1), nype_yz(0:nproc-1)
 integer, intent(out) :: nzpi_yz(0:nproc-1), nzpe_yz(0:nproc-1)
 integer, intent(out) :: nxpi_xz(0:nproc-1), nxpe_xz(0:nproc-1)
 integer, intent(out) :: nypi_xz(0:nproc-1), nype_xz(0:nproc-1)
 integer, intent(out) :: nzpi_xz(0:nproc-1), nzpe_xz(0:nproc-1)
 integer, intent(out) :: ksize(0:nproc-1)
 integer, intent(out) :: isize(0:nproc-1)
 integer, intent(out) :: jsize(0:nproc-1)
 integer, intent(out) :: nxpi_x(0:nproc-1), nxpe_x(0:nproc-1)
 integer, intent(out) :: nypi_y(0:nproc-1), nype_y(0:nproc-1)

 integer :: nstx, nsty, nstz
 integer :: rrank
 integer :: ierr
 integer :: coordn(1:3)
 integer :: i, j, k
 integer :: surplus

 if (decomp == 1) then

   nstx = int(float(nx)/float(dims(1)))
   nsty = int(float(ny)/float(dims(2)))
   nstz = int(float(nz+2)/float(dims(1))) ! Assumed dims(1) = dims(2)

   do j = 0, dims(2) - 1
     do i = 0, dims(1) - 1
       coordn(1) = i
       coordn(2) = j
       coordn(3) = 0
       call MPI_CART_RANK(COMM3D,coordn,rrank,ierr)
       nxpi(rrank) = coordn(1) * nstx + 1
       nxpe(rrank) = (coordn(1) + 1) * nstx
       nypi(rrank) = coordn(2) * nsty + 1
       nype(rrank) = (coordn(2) + 1) * nsty
       if (coordn(1) == 0)           nxpi(rrank) = 0
       if (coordn(1) == dims(1) - 1) nxpe(rrank) = nx + 1
       if (coordn(2) == 0)           nypi(rrank) = 0
       if (coordn(2) == dims(2) - 1) nype(rrank) = ny + 1 
       nzpi(rrank) = 0
       nzpe(rrank) = nz + 1
     enddo
   enddo
   
   nstx = int(float(nx)/float(dims(1)*dims(2)))
   nsty = int(float(ny)/float(dims(1)*dims(2)))
   do rrank = 0, nproc-1
     nxpi_x(rrank) = rrank*nstx + 1
     if (rrank == 0) nxpi_x(rrank) = 0
     nxpe_x(rrank) = (rrank + 1)*nstx
     if (rrank == nproc-1) nxpe_x(rrank) = nx + 1
     nypi_y(rrank) = rrank*nsty + 1
     if (rrank == 0) nypi_y(rrank) = 0
     nype_y(rrank) = (rrank + 1)*nsty
     if (rrank == nproc-1) nype_y(rrank) = ny + 1   
   enddo

 elseif (decomp == 2) then

   nstx = int(float(nx)/float(dims(1))) + 1
   nsty = int(float(ny)/float(dims(2))) + 1
   nstz = int(float(nz+2)/float(dims(1))) ! Assumed dims(1) = dims(2)

   do j = 0, dims(2) - 1
     do i = 0, dims(1) - 1

       coordn(1) = i
       coordn(2) = j
       coordn(3) = 0
       call MPI_CART_RANK(COMM3D,coordn,rrank,ierr)
  
       surplus = nstx*dims(1) - nx
       nxpi(rrank) = 0
       nxpe(rrank) = nstx - int(float(surplus)/2.)
       if (coordn(1) > 0) then
         do k = 1, coordn(1)
           nxpi(rrank) = nxpe(rrank) + 1
           nxpe(rrank) = nxpe(rrank) + nstx
         enddo
       endif
       if (coordn(1) == dims(1) - 1) nxpe(rrank) = nx + 1

       surplus = nsty*dims(2) - ny
       nypi(rrank) = 0
       nype(rrank) = nsty - int(float(surplus)/2.)
       if (coordn(2) > 0) then
         do k = 1, coordn(2)
           nypi(rrank) = nype(rrank) + 1
           nype(rrank) = nype(rrank) + nsty
         enddo
       endif
       if (coordn(2) == dims(2) - 1) nype(rrank) = ny + 1

       nzpi(rrank) = 0
       nzpe(rrank) = nz + 1

     enddo
   enddo

 endif


 do i = 0, dims(1) - 2
   ksize(i) =  nstz
 enddo
 ksize(dims(1)-1) = nz + 2 - (dims(1)-1)*nstz

 do i = 0, dims(1) - 1
   coordn(1) = i
   coordn(2) = 0
   call MPI_CART_RANK(COMM3D,coordn,rrank,ierr)
   isize(i) = nxpe(rrank) - nxpi(rrank) + 1
 enddo

 do i = 0, dims(2) - 1
   coordn(1) = 0
   coordn(2) = i
   call MPI_CART_RANK(COMM3D,coordn,rrank,ierr)
   jsize(i) = nype(rrank) - nypi(rrank) + 1
 enddo

 do j = 0, dims(2) - 1
   k = 0
   do i = 0, dims(1) - 1
     coordn(1) = i
     coordn(2) = j
     coordn(3) = 0
     call MPI_CART_RANK(COMM3D,coordn,rrank,ierr)
     nxpi_yz(rrank) = 0
     nxpe_yz(rrank) = nx + 1
     nypi_yz(rrank) = nypi(rrank)
     nype_yz(rrank) = nype(rrank)
     nzpi_yz(rrank) = k
     nzpe_yz(rrank) = k + ksize(i) - 1
     k = k + ksize(i)
   enddo
 enddo

 do i = 0, dims(1) - 1
   k = 0
   do j = 0, dims(2) - 1
     coordn(1) = i
     coordn(2) = j
     coordn(3) = 0
     call MPI_CART_RANK(COMM3D,coordn,rrank,ierr)
     nxpi_xz(rrank) = nxpi(rrank)
     nxpe_xz(rrank) = nxpe(rrank)
     nypi_xz(rrank) = 0
     nype_xz(rrank) = ny + 1
     nzpi_xz(rrank) = k
     nzpe_xz(rrank) = k + ksize(j) - 1
     k = k + ksize(j)
   enddo
 enddo

 END SUBROUTINE POIS_BOUNDS

